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Note: This is a working paper which will be expanded/updated frequently. One way to 
think about this paper is as an update of De Leeuw (2014), using simpler and more direct 
majorization methods. The directory deleeuwpdx.net/pubfolders/rstress has a pdf copy of 
this article, the complete Rmd hie that includes all code chunks, and R hies with the code. 

1 Problem 

Define the multidimensional scaling (MDS) loss function 

n 

a r (x) = ~ ( x'Aix ) r ) 2 , ( 1 ) 

i =1 

with r > 0 and the Aj positive semi-dehnite. The Wi are positive weights , the Si are non¬ 
negative dissimilarities. We call this rStress. Special cases are stress (Kruskal 1964) for 
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r = |, sstress (Takane, Young, and De Leeuw 1977) for r = 1, and the loss function used in 
MULTISCALE (Ramsay 1977) for r -y 0. 

Usually MDS is formulated uses Euclidean distances djfiX) between points j and £, which are 
rows of an n x p configuration matrix X. This fits into our formulation by setting x = vec(X) 
and by setting the to rip x np matrices of the form I p (g) E j( >. where the matrix E jf - has 
elements (j,j) and (£,£) equal to +1 and elements (j,£) and (£,j) equal to —1. Thus the 
are direct sums of p copies of the Ejg. Now x'AiX = dje(X). 

The problem we are trying to solve is to find an algorithm to minimize ay over x in for 
all values of r > 0. It is understood that parts of the algorithm may be different for different 
values of r. 


2 Minorization and Majorization 

We will design a convergent iterative majorization algorithm. Briefly this means constructing 
for each r > 0 a majorization function 7 r on M m ® M m such that a fix) < ^fix^y) for all x 
and y and such that a fix) = 7 fix,x) for all x. One step of the iterative algorithm is 

^(fc+i) _ argnun7 r (x, x®), 

unless xfi k ) already minimizes 7 r (x,x^), in which case we stop. If we do not stop we have 
the sandwich inequality 

afix {k+l) ) < 7 r (x (fc+1 ) ,x (fc) ) < lr [x {k \x {k) ) = afix {k) ). 

Thus majorization algorithms exhibit monotone convergence of loss function values. 

The history of majorization algorithms is complicated. They were first developed in the 
specific contexts of step-size estimation in descent algorothms (Ortega and Rheinboldt 1970), 
maximum likelihood estimation with missing data (Dempster, Laird, and Rubin 1977), and 
multidimensional scaling (De Leeuw 1977). Subsequently they were presented as a general 
class of optimization methods, and as a special case of block relaxation and augmentation 
methods, by De Leeuw (1994), see also Heiser (1995). The material in De Leeuw (1994) is 
(slowly, slowly) expanded into an e-book, with one part on block relaxation, augmentation 
and alternating least squares (“Homogeneity Analysis” 2015), one part on majorization (De 
Leeuw 2015a), and one part on mathematical background (De Leeuw 2015b). There have been 
many applicatons of the general majorization principle by The Dutch School in psychometrics, 
which includes Ten Berge, Kiers, Heiser, Mculman, Groenen, and Vermunt. 

Systematic use of majorization in statistics started with Lange, Hunter, and Yang (2000). 
There have been many further reviews, developments, and applications by Ken Lange and his 
co-workers. Theyo use the name MM Algorithms , where MM stands for either majorization- 
minimization or minorization-maximization. Lange (2013) has a very nice chapter on the 
MM algorithm, and he is currently working on an MM methods book. 
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3 Use of Homogeneity 

We start by exploring the obvious fact that 

min a r (x) = min min a(9x). 

X v ' 0>o x'x=l 

Let us introduce some notation to simplify the discussion. Without loss of generality we 
assume the dissimilarities are scaled by 


Y. w * S i = L 

2=1 


Next, it is convenient to define 


n 

Pr(x ) := J2 W i S i( X 'Ai x Y, 
2=1 


and 

n 

Vr( x ) ’■= J2 W i(x'AiX) 2r , 

2=1 

so that 


Now 


<J r (x) = 1 — 2 p r (x) + T]r(x). 


( 2 ) 


a r (9x) = 1 — 2 6 2r p r (x) + 9 4r r] r (x). 

Thus, as in De Leeuw (1977), thinking of rStress as a function of both x and 9 , 


min a r (9,x) = 1 — \ , 

o 9r(xy 

where the minimium is attained at 


( 3 ) 

( 4 ) 


2r = p r {;x) 

Vr(x) ' 


( 5 ) 


Thus minimizing cr r can be done by maximizing the ratio in (4) over x'x = 1. This is used 
in the majorization method proposed by De Leeuw (2014), by combining Dinkelbach (1967) 
and quadratic majorization. 

The approach of De Leeuw (2014) is somewhat cumbersome, because we first use homogeneity 
to reduce the MDS problem to a fractional programming problem, and then we use Dinkelbach 
to get rid of the fractional objective funcion again. Moreover one of the special cases requires 
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us to solve a modified eigenvalue problem of the type discussed, for example, by Hager (2001) 
in each iteration. 

In this paper we combine majorization with alternating least squares, working directly with 
(3). Thus we think of rStress as a function of both 9 and x, and we alternate minimization 
over 6 for fixed x and over x for fixed 6. Thus 

9^ = argnuncr r ($, x^), 

x ( k + l ) — ar g m j n a r (9 ( ' k \x). 

X f X=l 

The first step is trivially solved using (5). Updating x, however, is far from trivial and needs 
the majorization machinery we will develop in this paper. Since updating x is iterative, we 
have to choose how many majorization steps to take for updating x, before going to the next 
update for 9. 

Thus the main focus of the paper will be a majorization method for minimizing 


<j r (x, a) := 1 — 2 ap r (x) + a 2 r} r (x) (6) 

over x'x = 1 with a := 9 2r fixed at its current value. It is clear from (6) that we can hnd a 
majorization of a r by combining a minorization of p r and a majorization of r] r . 

4 Powers of Quadratic Forms 

We start with some lemmas we will need to construct our minorizations and majorizations. 
Lemma 4.1: f r (x ) := (x'Ax) r is convex on x'Ax > 0 if and only if r > 

Proof: The first and second derivative are 

Vf r {x) = 2 r(x' AxY^Ax, 


and 

V 2 f r (x) = 2 r(x'Ax) r ~ 1 (a + 2(r — \ ■ 

\ X /\.X J 

The matrix H r (x) := A + 2(r — 1) is psd for r = and its eigenvalues increase with r. 
Thus it is psd for all r > |. 

Also, if 0 < r < \ then, by Sylvester’s Law of Inertia, T> 2 f r (x) has precisely one negative 
eigenvalue, as well as rank(/l) — 1 positive eigenvalues, and n — rank(/l) zero eigenvalues. 
Thus in this case f r is not convex (and not concave either). QED 

By the way, we can use the definition of convexity to show lemma 4.1 is true for all x, not 
just for x with x'Ax > 0. 

Lemma 4.2: If r > | then 

f r {x) > (1 - 2 r)f r (y) + 2rf r -i(y)y'Ax. 
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Proof: Follows from the fact that f r is convex, i.e. above each of its tangents. QED 

Now write A(X) or Ax for the largest eigenvalue of a matrix X , and A(X) or Ax for the 
smallest eigenvalue. Note that if A = / (8) Eje then A^ = 2. 

Lemma 4.3: If r > 1 then 

A {V 2 f r (x)) < 2r(2r — 1)A^ (x'x) r ~ l . 


If x'x < 1 then 

X(V 2 f r (x)) < 2r(2r - 1)A^. 

Proof: If r > 1, then 

(u'Ax) 2 

u'H r (x)u = v!Au + 2(r — 1)---— < (2r — 1 )u Au. 

x Ax 


Thus 


X(H r (x)) < ( 2 r - 1)X A , 

and 

A (V 2 f r (x)) < 2r(2r — ^Xa^x'AxY^ 1 < 2r(2r — 1)A^ (x'a;) r_1 . 
Finally, if x'x < 1 then ( x'x) r ~ l < 1. QED 

Lemma 4.4: If 0 < r < 1 then 


fr(x) < (1 - r)f r (y) + rf r _ 1 (y)x , Ax. 


Proof: If r < 1 then ( x'Ax) r is concave in x'Ax , although not in x. Thus f r is below its 
tangents, and 

fr(x) < fr(y ) + r(y , Ay) r ~ 1 (x , Ax - y'Ay ), 
which simplifies to the required result. QED 
Lemma 4.5: If 0 < r < 1 and x'x = y'y = 1 then 

fr(x) < (1 - 2 r)f r (y) + 2rf r _ 1 (y)(x'(A - A A I)y + X A ). 


Proof: From 

x'Ax = y'Ay + 2y'A(x — y) + (x — y)'A(x — y) < —y'Ay + 2y'Ax + X A {x — y)'(x — y) 
we see that if x'x = y'y = 1 that 

x'Ax < 2x'(A — X A I)y + 2A^ — y'Ay. 

Substitute this in lemma 4.4 and simplify. QED 
Lemma 4.6: If 0 < r < | then 

A (V 2 f r (x)) > 2r(2r — l)X!Y A {x'x) r ~ l . 
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If x'x < 1 then 


k(V 2 f r {x)) > 2r(2r - 1)A^. 


Proof: We have 


as before. Thus 


(u' Ax) 2 
x' Ax 


< u'Au 


u'H{x)u > (2 r — 1 )u'Au > (2 r — 1)A au'u. 

The result follows because in addition x'Ax < A ax'x, and consequently (V Ax) r ~ l > 
A^ l {x'x) r ~ l . Moreover if x'x < 1 then ( x'x ) r > 1. QED 

The following lemma, defining a type of uniform quadratic majorization (De Leeuw and 
Lange 2009), is an additional useful tool. Because we work on the unit sphere, the quadratic 
majorization actually becomes linear. This relies heavily on the fact that if x'x = y'y = 1 
and z is on the segment [x,y] connecting x and y, then z'z < 1. 

Lemma 4.7: Suppose <f> is any homogeneous function of degree s, x'x = y'y = 1 and 
A (V 2 <p(z)) < k for all z G [x,y]. Then 

4>{x) < (1 - s)(f(y) +n + x'{V<j>{y) - ny). 

In the same way, if A {V 2 (j)(z)) > k for all z G [x,y\ we have 

4>{x) > (1 - s)<f>{y) + k + x'(V<j)(y) - ny). 

Proof: We only prove the first part. The proof of the second part goes the same. By Taylor’s 
theorem we have for all x and y 


<t>(x) < <P(v) + (x ~ y)"D<p(y) + -k(x - y)'(x - y), 
which simplifies to the stated result if x'x = y'y = 1 and <f> is homogeneous of degree s. QED 


5 Majorizing rStress 

We will now use the inequalities from the lemmas in the previous section to derive linear 
minorizations of a r , by combining majorization of r/ r and minorization of p r . Our majorizations 
are in the form f(x) < x'g(y) + e(y), where both g and e are functions which do not depend 
on x , only on y. Since the precise form of e has no influence on the iterative algorithm we 
usually leave it unspecified, so in different formulas e can refer to different funcions. Clearly 
the majorization algorithm will have the form 


r (k+D _ y( x(k) ) 


For the two key results of this paper we need to define matrices 


( 7 ) 
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( 8 ) 


B r(y) ■= Yj w My'AivY lA i, 

i =1 

and 

n 

c r (y) ■=J2 w i(y' A iy) 2r l Ai, ( 9 ) 

%— l 

and scalars 

n 

Pr := (2r - 1) wA^( A i), (10) 

1=1 

and 

n 

i r ■= J2 w i(y' A iy) 2r ~ l K A i), (n) 

i= 1 

and 

n 

S r := (4r - 1) ^WjA 2r (A)- (12) 

i=l 

Note that S r = (3 2r - Also y'B r (y)y = p r (y) and y'C r (y)y = r/ r (y). 

Theorem 5.1: If 0 < r < \ then (7) with 

g(y) = ~{{B r (y ) - &/) - a(C r (y) - 7 V-Oh/ (13) 

is a convergent majorization algorithm. 

Proof: To majorize rj r we use lemma 4.5. After some calculation we find the linear majoriza¬ 
tion 

Vr(x) < Arx’(C r (y ) - 7 r /)y + e(y), 
where e is some function which only depends on y and not on x. 

To minorize p r we use lemma 4.7 with <p equal to f r . Note that f r is homogeneous of order 
2r. Thus 

(x'Aix) r > (1 - 2r)(y , A i y) r + m + x , (2r(y , A i y) r ~ 1 A i - K t I)y 
By lemma 4.6 we have K t = 2r(2r — l)A r (Aj). Thus 

p r (x) > 2 rx\B r (y) - p r I)y + e(y). 
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By combining the results for p r and g r we see that 

cr r (x) < 4 arx'giy) + e(y) 

with g(y ) dehned by (13). QED 
Theorem 5.2: If r > | then (7) with 

g(y) : = ~{( B r(y) - a{C r (y) - S r i)}y (14) 

is a convergent majorization algorithm. 

Proof: We can use lemma 4.2 to get 

p r (x) > 2 rx'B r (y)y + e(y). 

Because p r is homogeneous of order 4r we see from lemma 4.7 that 

p r (x) < 4 rx'(C r (y) - S r I)y + e(y). 

By combining the results for p r and r/ r we see that 

cr r (x) < Arax’g(y) + e(y) 

with g(y) dehned by (14). QED 


6 Special Considerations 


If' r — f we have /3 r = 0, and 


B ' {y)= tr‘vm 


A-i, 


and 

n 

Cr(y) = J2 W i A i- 

i =1 

These are precisely the matrices used in the SMACOF algorithm of De Leeuw (1977), 
implemented in De Leeuw and Mair (2009), which is 


x ( k +i) = c+(x w )B r (x {k) )x {k \ 


where C*+ is the Moore-Penrose inverse. Also note that if r = | then -y r = S r and the 
algorithms of theorems 5.1 and 5.2 are the same. 

The case r —>• 0 also requires some special attention. Define 


do Or) 


X>i(log<*i - log x'Aix) 2 . 

i —1 
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Note that 


53 lUi(log ^ - log {x'Aix) r f = r 2 a 0 (x), 

i=i 

which means that minimizing JYi=i w i(l°g Y ~ log (a; 7 Aix) r ) 2 for any r > 0 can be done by 
minimizing cr 0 . 

Now use 

x r — 1 

log x = lim-. 

r—» 0 f 


Thus, for small r, 


and 


log Si — log x'AiX 


S\ — (. x'Aix) r 


cr 0 {x) ~ ~^J2wi(Sl - (. x'AixYf 


i —1 


On the other hand, if S % « x'A^ then for any r we have the approximation 


log {x'Aix) T fa log SI + —((x'AixY - <$[), 


so that 

" w- 

M x ) ~ ~ ( x 'Ax) r ) 2 . 

i =1 

Both approximatations to oo can be minimized by our iterative majorization algorithm for 

r < - 

i ^ 2 - 

Note that our definition of rStress compares dissimilarities with powers of Euclidean distances. 
Alternatively, we can compare powers of dissimilarities with corresponding powers of Euclidean 
distances, i.e. we can define another rStress as 


cr r (x) := ^ Z w i((SiY ~ ( x'AixY ) 2 - (15) 

i =1 

This is similar to our approach in the limiting case r —> 0. We do not need an additional 
algorithm to minimize ay, because we can just input the (Sf) r instead of S t in our previous 
majorization method. 

The difference between ay and ay disappears, of course, if we make the algorithm nonmetric. In 
that case we alternate minimization over x with minimization over monotone transformations 
of S (or of 5 r , which is of course the same thing). We could easily add a monotone regression 
step to our alternating least squares algorithm. 
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7 Examples 

7.1 Artificial 

Let’s start with a small example that shows the flow of the algorithm. We perform one 
iteration and write out the intermediary results in detail. 

Suppose n = 3, the weights w are all one, and the dissimilarities are 1, 2, 3. The matrices A 
are 

## C, 1] [, 2] [, 3] [ ,4] 

## [ 1 ,] 1 1 - 1-1 

## [ 2 ,] 1 1 - 1-1 

## [3,] -1-1 1 1 

## [4,] -1-1 1 1 

## C, 1] [, 2] [, 3] [,4] 

## [1,] 1-1-1 1 

## [ 2 ,] -1 1 1-1 

## [3,] -1 1 1-1 

## [4,] 1-1-1 1 

## C, 1] C, 2] [,3] [,4] 

## [ 1 ,] 1-1 1-1 

## [ 2 ,] -1 1-1 1 

## [3,] 1-1 1-1 

## [4,] -1 1-1 1 

The largest eigenvalues of all three A matrices are equal to 4. Thus S r of (??) is 144. 

Let’s start with x four random normals. We want to minimize rStress for r = 1, i.e. sStress. 
## [1] 3.4514035 0.2191467 0.0485130 

At the start sStress is 17.892093. Also p is 4.0352358 and rj is 11.9625647. This means a is 
0.337322 and the minimum of rStress over a is 12.6388263. 

We now compute the matrices B r and C r at x. 

For B r we have 


## 

## 

[1J 

[,1] 

6 

[, 2] 
-4 

1-1 

GO 
O >— 1 

[ >4] 
-2 

## 

[2,] 

-4 

6 

-2 

0 

## 

[3,] 

0 

-2 

6 

-4 

## 

[4,] 

-2 

0 

-4 

6 


and for C r 

## C,1] [,2] [,3] [,4] 

## [1,] 3.719063 3.183744 -3.622037 -3.280770 
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## [2,] 3.183744 3.719063 -3.280770 -3.622037 

## [3,] -3.622037 -3.280770 3.719063 3.183744 
## [4,] -3.280770 -3.622037 3.183744 3.719063 

Also, S r is equal to 144 and thus 

Now g is 

## [1] 9.408613 11.604293 -1.162210 -7.853545 

and the new x is 

## [1] -0.55613830 -0.68592383 0.06869765 0.46421904 

With this x sStress is 12.4979159. We now go back, compute a new a, and so on. 


7.2 Dutch Political Parties 

De Gruijter (1967) collected dissimilarity measures for nine Dutch political parties. 

## KVP PvdA VVD ARP CHU CPN PSP BP 

## PvdA 5.63 

## VVD 5.27 6.72 

## ARP 4.60 5.64 5.46 

## CHU 4.80 6.22 4.97 3.20 

## CPN 7.54 5.12 8.13 7.84 7.80 

## PSP 6.73 4.59 7.55 6.73 7.08 4.08 

## BP 7.18 7.22 6.90 7.28 6.96 6.34 6.88 

## D66 6.17 5.47 4.67 6.13 6.04 7.42 6.36 7.36 

The solutions for 0.1, 0.25, 0.5, 0.75, 1, 2 are given in figure 1. 
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Figure 1: De Gruijter Data: Separate Configurations 

The basic structure for all solutions is the same, we see the familiar horseshoe , in which the 
left-right scale is bended so that extreme left (CPN) and extreme right (BP) are close. This 
is perhaps most clear for r = .75. If r becomes large, then we see increasing clustering. 

We also give all six solutions plotted on top of each other (after Procrustus matching) in 
figure 9. It shows that the two extreme parties are stable over solutions, but the middle is 
more variable. 
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dimension 1 

Figure 2: De Gruijter Data: Joint Configuration 

Plotting rStress in figure 3 as a function of r shows, somewhat surprisingly perhaps, an 
increasing function. Because the data are average dissimilarities, it is likely there is a fairly 
large additive constant, and additive constants correspond with lower values of r. 



r 

Figure 3: De Gruijter Data: rStress 

We put the nuber of iterations and the rStress values in a small table. 
## r: 0.10 iterations: 29103 rStress: 0.005464 
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## 

r: 

0.25 

iterations: 

3605 

rStress: 

0.006310 

## 

r: 

0.50 

iterations: 

3566 

rStress: 

0.044603 

## 

r: 

0.75 

iterations: 

3440 

rStress: 

0.107113 

## 

r: 

1.00 

iterations: 

100000 

rStress: 

0.155392 

## 

r: 

2.00 

iterations: 

100000 

rStress: 

0.234877 


Figure 4 shows the six Shepard diagrams. We see the clustering for r = 2. For r = .1 the 
Shepard diagram becomes concave, indicating that the larger dissimilarities are underestimated 
and reflecting the fact that for small powers the powered distances will all be close to one. 
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Figure 4: De Gruijter Data: Shepard Diagrams 


7.3 Ekman Color Data 

We use the color data from Ekman (1954), with two dimensions and unit weights. As we 
know from previous analyses, MDS of the Ekman data gives a very good fit and a very clear 
representation of the color circle. 

## 434 445 465 472 490 504 537 555 584 600 610 628 651 

## 445 0.14 

## 465 0.58 0.50 

## 472 0.58 0.56 0.19 

## 490 0.82 0.78 0.53 0.46 

## 504 0.94 0.91 0.83 0.75 0.39 

## 537 0.93 0.93 0.90 0.90 0.69 0.38 
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## 555 0.96 0.93 0.92 0.91 0.74 0.55 0.27 

## 584 0.98 0.98 0.98 0.98 0.93 0.86 0.78 0.67 

## 600 0.93 0.96 0.99 0.99 0.98 0.92 0.86 0.81 0.42 

## 610 0.91 0.93 0.98 1.00 0.98 0.98 0.95 0.96 0.63 0.26 

## 628 0.88 0.89 0.99 0.99 0.99 0.98 0.98 0.97 0.73 0.50 0.24 

## 651 0.87 0.87 0.95 0.98 0.98 0.98 0.98 0.98 0.80 0.59 0.38 0.15 

## 674 0.84 0.86 0.97 0.96 1.00 0.99 1.00 0.98 0.77 0.72 0.45 0.32 0.24 

In our algorithm we always start with a classical metric scaling solution (Torgerson 1958). 
We minimize rStress for r equal to 0.1, 0.25, 0.5, 0.75, 1, 2. The six solutions are plotted 
jointly in Figure 5. Because of the very good fit, solutions for different values of r are similar 
and very much on the color circle. 



Figure 5: Ekman Color Data: Configuration 

The Shepard diagrams in figure 6 are interesting, because, unlike the configurations, they are 
quite different, becoming more and more convex as r increases. The plot for r = .25 looks 
best. 
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Figure 6: Ekman Color Data: Shepard Diagrams 


We plot rStress as a function of r in Figure 7. The best fit is attained for r = .25, which 
means fitting square roots of Euclidean distances to dissimilarities. Nonmetric MDS of the 
Ekman data has already indicated that the optimal nonmetric fit is attained with a concave 
increasing transformation. 
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Figure 7: Ekman Color Data: rStress 


## 

r: 

0.10 

iterations: 

100000 

rStress: 

0.017839 

## 

r: 

0.25 

iterations: 

1361 

rStress: 

0.001910 

## 

r: 

0.50 

iterations: 

535 

rStress: 

0.017213 

## 

r: 

0.75 

iterations: 

3343 

rStress: 

0.054769 

## 

r: 

1.00 

iterations: 

13749 

rStress: 

0.093063 

## 

r: 

2.00 

iterations: 

100000 

rStress: 

0.181719 


Especially when r is far from | it will make a difference if we minimize oy or the ay of (15). 
We illustrate this by choosing r = .01, which will presumably take us close to the logarithm. 
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dimension 1 

Figure 8: Ekman Color Data: Power Transformation 

After 14837 iterations we find a ay value of 0.000012. We see a rather clear deviation from 
the circular pattern we found with rStress. 


8 Code 

torgerson <- function(delta, p = 2) { 
doubleCenter <- function(x) { 
n <- dim(x) [1] 
m <- dim(x) [2] 
s <- sum(x) / (n * m) 
xr <- rowSums(x) / m 
xc <- colSums(x) / n 
return((x - outerfxr, xc, "+")) + s) 
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} 

z <- 

eigen(-doubleCenter ( (as.matrix (delta) 2) / 2) , symmetric = TRUE) 
v <- pmax(z$values, 0 ) 

return(z$vectors[, 1 : p] %*% diag(sqrt (v[ 1 :p]))) 

} 

enorm <- function (x, w = 1 ) { 

return (sqrt (sum (w * (x ~ 2 )))) 

} 

sqdist <- function (x) { 
s <- tcrossprod (x) 
v <- diag (s) 

return (outer (v, v, "+") - 2 * s) 

} 

mkBmat <- function (x) { 
d <- rowSums (x) 
x <- -x 
diag (x) <- d 
return (x) 


mkPower <- function (x, r) { 
n <- nrow (x) 

return (abs ((x + diag (n)) r) - diag(n)) 

} 

matchConf <- function (x, 

eps = le-6, 
itmax = 100, 
verbose = TRUE) { 

m <- length (x) 
n <- nrow (x[[1]]) 
p <- ncol (x[[1]]) 
itel <- 1 
dold <- Inf 
repeat { 

y <- matrix (0, n, p) 
for (k in l:m) 
y <- y + x [ [k] ] 
y <- y / m 
dnew <- 0 
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for (k in l:m) { 

s <- svd (crossprod(x [[k]], y)) 

x[[k]] <- tcrossprod (x[[k]j %*% (s$u) , s$v) 

dnew <- dnew + sum ((y - x[[k]]) ~ 2 ) 

> 

if (verbose) { 
cat ( 

formate (itel, width = 4, format = "d") , 
formate ( 
dold, 

digits = 10, 
width = 13, 
format = "f" 

), 

formate ( 
dnew, 

digits = 10, 
width = 13, 
format = "f" 

), 

"\n" 

) 

> 

if ((itel == itmax) I I ((dold - dnew) < eps)) 
break () 

itel <- itel + 1 
dold <- dnew 

} 

return (x) 


rStressMin <- 

function (delta, 

w = 1 - diag (nrow (delta)), 

p = 2, 

r = 0.5, 
eps = le-10, 
itmax = 100000, 
verbose = TRUE) { 
delta <- delta / enorm (delta, w) 
itel <- 1 

xold <- torgerson (delta, p = p) 
xold <- xold / enorm (xold) 
n <- nrow (xold) 
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nn <- diag (n) 
dold <- sqdist (xold) 

rold <- sum (w * delta * mkPower (dold, r)) 
nold <- sum (w * mkPower (dold, 2 * r)) 
aold <- rold / nold 

sold <- 1 - 2 * aold * rold + (aold 2 ) * nold 
repeat { 

pi <- mkPower (dold, r - 1) 
p2 <- mkPower (dold, (2 * r) - 1) 
by <- mkBmat (w * delta * pi) 
cy <- mkBmat (w * p2) 
ga <- 2 * sum (w * p2) 

be <- (2 * r - 1 ) * (2 ~ r) * sum (w * delta) 

de <- (4 * r - 1) * (4 ~ r) * sum (w) 

if (r >= 0.5) { 

my <- by - aold * (cy - de * nn) 

} 

if (r < 0.5) { 

my <- (by - be * nn) - aold * (cy - ga * nn) 

} 

xnew <- my °/ 0 *°/ 0 xold 

xnew <- xnew / enorm (xnew) 

dnew <- sqdist (xnew) 

rnew <- sum (w * delta * mkPower (dnew, r)) 
nnew <- sum (w * mkPower (dnew, 2 * r)) 
anew <- rnew / nnew 

snew <- 1 - 2 * anew * rnew + (anew 2 ) * nnew 
if (verbose) { 
cat ( 

formate (itel, width = 4, format = "d"), 
formate ( 
sold, 

digits = 10 , 
width = 13, 
format = "f" 

), 

formate ( 
snew, 

digits = 10 , 
width = 13, 
format = "f" 

), 

"\n" 

) 
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} 

if ((itel == itmax) II ((sold - snew) < eps)) 


break () 

itel 

<- itel + 1 

xold 

<- xnew 

dold 

<- dnew 

sold 

<- snew 

aold 

<- anew 

> 

return 

(list (x = xnew. 


alpha = anew, 
sigma = snew, 
itel = itel)) 


9 NEWS 


001 01/12/16 
002 01/12/16 
003 01/13/16 
004 01/13/16 
005 01/13/16 
006 01/13/16 
007 01/13/16 
008 01/13/16 
009 01/13/16 
010 01/14/16 
011 01/14/16 


First published version, without code and examples 
Added artificial example and code. 

Squashed two nasties 
Added Ekman Example 
Many small edits 
Reorganized proofs 
Added political parties example 
Added configuration matching 
Added Ekman power solution 
Text in examples - Additional plots 
small change in theorems and code 
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